﻿using System;
using System.IO;

namespace NumericalSolution {
  class Diffusion1DFEMMain {

    static String outfile = "c:/nsworkspace/diffusion1dfem.txt";

    public static void Main() {
      Diffusion1DFEM sim = new Diffusion1DFEM();
      sim.NumberOfElements = 21;
      sim.NumberOfNodesInAnElement = 2;
      sim.DeltaX = 1.0 / sim.NumberOfElements;
      sim.NumericalIntegration = new GaussLagrange(2);
      sim.EndTime = 1.0;
      sim.DeltaTime = 0.005;
      sim.CreateMesh();
      sim.Solve();

      double[] u;
      StreamWriter writer = File.CreateText(outfile);

      double error = 0.0;
      double count = 0;
      double maxError = 0.0;
      for (sim.CurrentTime = 0.0;
           sim.CurrentTime <= sim.EndTime;
           sim.CurrentTime += sim.DeltaTime) {

        if (sim.CurrentTime == 0.0) {
          u = sim.Initialize();
        } else {
          u = sim.Next();
        }

        for (int i = 0; i < u.Length; i++) {
          double t = sim.CurrentTime;
          double x = sim.DeltaX * i;
          double n = u[i];
          double e = exact(x, t);
          writer.Write(t + " ");
          writer.Write(x + " ");
          writer.WriteLine(n);

          if (e != 0) {
            error = error + Math.Abs(n - e) / e;
            count++;
            double te = Math.Abs(n - e);
            if (te > maxError) {
              maxError = te;
            }
          }
        }
        writer.WriteLine();
      }
      error /= count;
      writer.WriteLine("%ERROR: " + error);
      writer.WriteLine("%MAX ERROR: " + maxError);
      writer.Close();
    }

    public static double exact(double x, double t) {
      return Math.Sin(Math.PI * x) * Math.Exp(- Math.PI * Math.PI * t);
    }
  }
}
